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Numerical exploration of the properties of singularities could, in principle, yield 
detailed understanding of their nature in physically realistic cases. Examples of nu- 
merical investigations into the formation of naked singularities, critical behavior in 
collapse, passage through the Cauchy horizon, chaos of the Mixmastcr singularity, 
and singularities in spatially inhomogeneous cosmololgies are discussed. 



1 Introduction 

The singularity theorems BiSi state that Einstein's equations will not evolve 
regular initial data arbitrarily far into the future or the past. An obstruc- 
tion such as infinite curvature or the termination of geodesies will always 
arise to stop the evolution somewhere. The simplest, physically relevant solu- 
tions representing for example a homogeneous, isotropic universe (Friedmann- 
Robertson- Walker (FRW)) or a spherically symmetric black hole (Schwarz- 
schild) contain space-like infinite curvature singularities. Although, in princi- 
ple, the presence of a singularity could lead to unpredictable measurements for 
a physically realistic observer, this does not happen for these two solutions. 
The surface of last scattering of the cosmic microwave background in the cos- 
mological case and the event horizon in the black hole (BH) case effectively 
hide the singularity from present day, external observers. The extent to which 
this "hidden" singularity is generic and the types of singularities that appear 
in generic spacetimes remain major open questions in general relativity. The 
questions arise quickly since other exact solutions to Einstein's equations have 
singularities which are quite different from those described above. For exam- 
ple, the charged BH (Reissner-Nordstrom solution) has a time-like singularity. 
It also contains a Cauchy horizon (CH) marking the boundary of predictability 
of space-like initial data. A test observer can pass through the CH to another 
region of the extended spacetime. More general cosmologies can exhibit singu- 
larity behavior different from that in FRW. The Big Bang in FRW is classified 
as an asymptotically velocity term dominated (AVTD) singularity Era since any 
spatial curvature term in the Hamiltonian constraint becomes negligible com- 
pared to the square of the expansion rate as the singularity is approached. 
Howeyfijj, some anisotropic, homogeneous models exhibit Mixmaster dynamics 
(MD)Qla and are not AVTD — the influence of the spatial scalar curvature can 
never be neglected. 

Once the simplest, exactly solvable models are left behind, understanding 
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of the singularity becomes more difficult. Other chapters BO describe recent 
analytic progress. However, such methods yield either detailed knowledge of 
unrealistic, simplified (usually by symmetries) spacetimes or powerful, general 
results that do not contain details. To overcome these limitations, one might 
consider numerical methods to evolve realistic spacetimes to the point where 
the properties of the singularity may be identified. Of course, most of the effort 
in numerical relativity applied to BH collisions has addressed the avoidance 
of singularities El One wishes to keep the computational grid in the observable 
region outside the horizon. Much less computational effort has focused on the 
nature of the singularity itself. Numerical calculations, even more than analytic 
ones, require finite values for all quantities. Ideally then, one must describe 
the singularity by the asymptotic non-singular approach to it. A numerical 
method which can follow the evolution into this asymptotic regime will then 
yield information about the singularity. Since the numerical study must begin 
with a particular set of initial data, the results can never have the force of 
mathematical proof. One may hope, however, that such studies will provide 
an understanding of the "phenomenology" of singularities that will eventually 
guide and motivate rigorous results. 

In the following, we shall consider examples of numerical study of singular- 
ities both for asymptotically flat (AF) spacetimes and for cosmological mod- 
els. These examples have been chosen to illustrate primarily numerical studies 
whose focus is the nature of the singularity itself. In the AF context, we shall 
consider two questions. The first is whether or not naked singularities exist 
for realistic matter sources. One approach is to explore highly non-spherical 
collapse looking for spindle or pancake singularities. If the formation of an 
event horizon requires a limit on the aspect ratio of the matterpl such configu- 
rations may yield a naked singularity. Another approach is to probe the limits 
between initial configurations which lead to black holes and those which yield 
no singularity at all (i.e. flat spacetime plus radiation) to explore the singular- 
ity as the BH mass goes to zero. This quest led-jiaturally to the discovery of 
critical behavior in the collapse of a scalar fieldEj The other question which is 
now beginning to yield to numerical attack involves the stability of the Cauchx 
horizon in charged or rotating black holes where it has been conjectured BEJ 
that a real observer, as opposed to a test mass, cannot pass through the CH 
since realistic perturbed spacetimes will convert the CH to a true singularity. 
In cosmology, we shall consider both the behavior of the Mixmaster model and 
the issue of whether or not its properties are applicable to generic cosmolog- 
ical singularities. Although numerical evolution of the Mixmaster equations 
has a long history, recent developments are motivated by inconsistencies be- 
tween the known sensitivity to initial conditions and standard measures of the 
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chaos usually associated with such behavior00B@fl@li] Bclinskii, Khalat- 
nikov, and Lifshitz (BKL) long ago claimed that it is possible to formulate 
the generic cosmological solution to Einstein's equations near the singularity 
as a Mixmaster universe at every point. While others have questioned the the 
validity of this claimpj there is as yet no evidence either way for such behavior 
in spatially inhomogeneous cosmologies. We shall discuss a numerical program 
to address this issue. 

2 Singularities in AF spacetimes 

2.1 Naked singularities and the hoop conjecture 

The strong cosmic censorship conjecture §1 requires a singularity formed from 
regular, asymptotically flat initial data to be hidden from an external observer 
by an event horizon. Counter examples have been known for a long time 
but tend to be dismissed as unrealistic in some way. The goal of a numeri- 
cal approach is to try to search for naked singularities arising from physically 
reasonable initial conditions. Ajaossible regime for such systems is motivated 
by Thome's "hoop conjecture" EJ that collapse will yield a black hole only if 
a mass M is compressed to a region with circumference C < 4nM in all di- 
rections. (Note that one must take care to define C and M especially if the 
initial data are not at least axially symmetric.) If the hoop conjecture is true, 
naked singularities may form if collapse can yield C > AirM in some direction. 
The existence of a naked singularity is inferred from the absence of an appar- 
ent horizon (AH) which can be identified locally by following null geodesies. 
Although a definitive identification of a naked singularity requires the event 
horizon (EH) to be proven to be absent, to identify an EH requires knowledge 
of the entire spacetime. Methods to find an EH in a numerically determined 
spacetime—have only recently become available and have not been applied to 
this iss.ueE3 To attempt to produce naked singularities, Shapiro and Teukolsky 
(ST) Ell considered collapse of prolate spheroids of collisionless gas. (Nakamura 
and Satot3 had previously studied the collapse of non- rotating deformed stars 
with an initial large reduction of internal energy and apparently found spindle 
or pancake singularities in extreme cases.) ST solved the general relativistic 
Vlasov equation for the particles along with Einstein's equations for the grav- 
itational field. Null geodesies were followed to identify an AH if present. The 
curvature invariant / = R^ t ,p <T R lluf " 7 was also computed. They found that an 
AH (and presumably a BH) formed if C < AttM < 1 everywhere but no AH 
(and presumably a naked singularity) in the opposite case. In the latter case, 
the evolution (not surprisingly) could not proceedjaast the moment of forma- 
tion of the singularity. In a subsequent study, STE3 also showed that a small 
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amount of rotation (counter rotating particles with no net angular momen- 
tum) does not jjrevent the formation of a naked spindle singularity. However, 
Wald and IyerEZl have shown that the Schwarzschild solution has a time slicing 
whose evolution approaches arbitrarily close to the singularity with no AH in 
any slice (but, of course, with an EH in the spacetime). This may mean that 
there is a chance that the increasing prolateness found by ST in effect changes 
the slicing to one with no apparent horizon just at the point required by the 
hoop conjecture. While, on the face of it, this seems unlikely, Tod gives an 
example where the AH does not form on a chosen constant time slice — but 
rather different portions form at different times. He argues that a numerical 
simulation might be forced by the singularity to end before the formation of 
the AH is complete. Such an AH would not be found by the simulationsES In 
response, Shapiro and Teukolsky considered equilibrium sequences of prolate 
relativistic star clustersE2l The idea is to counter the possibility that an EH 
might form after the simulation must stop. If an equilibrium configuration is 
non-singular, it cannot contain an EH since singularity theorems say that an 
EH implies a singularity. However, a sequence of non-singular equilibria with 
rising I ever closer to the spindle singularity would lend support to the exis- 
tence of a naked spindle singularity since one can approach the singular state 
without formation of an EH. They constructed this sequence and found that 
the singular end points were very similar to their dynamical spindle singularity. 

Motivated by ST's results^ EcheverriaEl numerically studied the proper- 
ties of the naked singularity that is known to form in the collapse of an infinite, 
cylindrical dust shelllij While the asymptotic state can be found analytically, 
the approach to it must be followed numerically. The analytic asymptotic solu- 
tion can be matched to the numerical one (which cannot be followed all the way 
to the collapse) to show that the singularity is strong (an observer experiences 
infinite stretching parallel to the symmetry axis and squeezing perpendicular 
to the symmetry axis). A burst of gravitational radiation emitted just prior 
to the formation of the singularity can be studied. 

A_yery recent numerical study of the hoop conjecture was made by Chiba 
ct alEEl Rather than a dynamical collapse model, they searched for AH's in 
analytic initial data for discs, annuli, and rings. Previous studies of this type 
were done bv_.Nakamura et al Eij with oblate and prolate spheroids and by 
WojtkiewiczLj with axisymmetric singular lines and rings. The summary of 
their results islhat an AH foxms if C < AirM < 1.26. (Analytic results due to 
Barrabes et al r 5 !! 36 ! and Todt3 give similar quantitative results with different 
initial data classes and (possibly) definition of C.) The results of all these 
searches for naked singularities are controversial but could be resolved if the 
presence or absence of the EH could be determined. 
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2.2 Critical behavior in collapse 

We now consider an effect that was originally found by ChoptuikEl in a numer- 
ical study of the collapse of a spherically symmetric massless scalar field. We 
note that this is the first completely new phenomenon in general relativity to 
be discovered by numerical simulation. In collapse of a scalar field, essentially 
two things can happen: either a BH forms or the scalar waves pass through 
each other and disperse. Choptuik discovered that for any 1-parameter set of 
initial data labeled by p, there is a critical value p* such that p > p* yields a 
BH. He found 

M B H~c F (p~ P y (i) 

where Mbh is the mass of the eventual BH. The constant Cf depends on 
the parameter of the initial data that is selected but 7 rj .37 is the same for 
all choices. Furthermore, in terms of logarithmic variables p — lnr + re, r = 
ln(Tp — To) + re (To is the proper time of an observer at r = with T * the finite 
proper time at which the critical evolution concludes), the waveform X repeats 
(echoes) at intervals of A in r if p is rescaled to p — A, i.e. X(p — A, r — A) rj 
X(p,r). The scaling behavior ([!]) demonstrates that the minimum BH mass 
(for bosons) is zero. The critical solution itself is a counter-example to cosmic 
censorship (since the formation of the zero mass BH causes high curvature 
regions become visible at r = 00). (See, e.g., the discussion in Hirschmann and 
EardleyB) 

Soon after this discovery, scaling andjnitical phenomena were found in a 
variety of contexts. Abrahams and EvanstS discovered the same phenomenon 
in axisymmetric gravitational wave collapse with a different value of A and, to 
within numerical error, the same value of 7. (Note that the rescaling of r with 
e A w 30 required Choptuik to use adaptive mesh refinement (AMR) to distin- 
guish subsequent echoes. Abrahams and Evans' smaller A (e A r> 1.8) allowed 
thep. to see echoing with their 2 + 1 code without AMR.) Recently, Garfin- 
kleE3 has confirmed Choptuik's results with a completely different algorithm 
that does not require AMB_ His use of Goldwirth and Piran's a method of 
simulating Christodoulou'sEH formulation of the spherically symmetric scalar 
field in null coordinates allowed the grid to be automatically rescaled by choos- 
ing the edge of the grid to be the null ray that just hits the central observer 
at the end of the critical evolution. (Missing points of null rays that cross the 
central observer's world line arejceplaced by interpolation between those that 
remain.) Hamade and Stewart Il3 have also repeated Choptuik's calculation 
using null coordinates and AMR. They are able to achieve greater accuracy 
and find 7 = .374. j—. 

Evans and Coleman c2l realized that self-similar rather than self-periodic 
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collapse might be more tractable both numerically (since ODE's rather than 
PDE's are involved) and analytically. They discovered that a collapsing radi- 
ation fluid had that desirable property. (Note that self-similarity (homothetic 
motion) is incompatible with AF. However, most of the action occurs in the 
center so that a match of the self-similar inner region to an outer AFj^n^should 
always be possible.) In a series of papers, Hirschmann and EardleyElI'Lj devel- 
oped a (numerical) self-similar solution to the spherically symmetric complex 
scalar field equations. These are ODE's with too many boundary conditions 
causing a solution to exist only for certain fixed values of A. Numerical so- 
lution of this eigenvalue problem allows very accurate determination of A. 
The self-similarity also allows accurate calculation of 7 as follows: The critical 
p = p* solution is unstable to a small change in p. At any time t (where t < 
is increasing toward zero), the amplitude a of the perturbation exhibits power 
law growth: 

a X (p-p*)(-t)- K (2) 

where k > 0. At any fixed t, larger a implies larger Mbh- Equivalcntly, 
any fixed amplitude a — S will be reached faster for larger eventual Mbh- 
Scaling arguments give the dependence of Mbh on the time at which any 
fixed amplitude is reached: 

M BH oc (-h) (3) 

where 

(p-p*)(-hT*5. (4) 

Thus 

Mmk(p-P*) 1/k . (5) 

Therefore, one need only identify the growth rate of the unstable mode to 
obtain an accurate value of 7 = 1/k. It is not necessary to undertake the 
entire dynamical evolution or probe the space of initial data. Hirschmann and 
Eardkiv obtain 7 = 0.387106 for the complex scalar field solution while Koiki 
et al C3 obtain 7 — 0.35580192 for the Evans-Coleman solution. Although 
the similarities among the critical exponents 7 in the collapse computations 
suggested a universal value, MaisoncEl used these same scaling-perturbation 
methods to show that 7 depends on the equation of state_p = kp of the fluid 
in the Evans-Coleman solution. Very recently, Gundlachcll used a similar ap- 
proach to locate Choptuik's critical solution accurately. This is much harder 
due to its discrete self-similarity. He reformulates the model as nonlinear hy- 
perbolic boundary value problem with eigenvalue A and finds A = 3.4439. As 
with the self-similar solutions described above, the critical solution is found di- 
rectly without the need to perform a dynamical evolution or explore the space 
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of initial data. Perturbation of this solution could yield an accurate value of 
7- 

2.3 Nature of the singularity in charged or rotating black holes 

Unlike the simple singularity structure of the Schwarzschild solution, where the 
event horizon encloses a spacelike singularity at r = 0, charged and/or rotating 
BH's have a much richer singularity structure. The extended spacetimes have 
an inner Cauchy horizon (CH) which is the boundary of predictability.. To 
the future of the CH lies a timelike (ring) singularity!]! Poisson and Israel Hc3 
began an analytic study of the effect of perturbations on the CH. Their goal 
was to check conjectures that the blue-shifted infalling radiation during collapse 
would convert the CH into a true singularity and thus prevent an observer's 
passage into the rest of the extended regions. By including both ingoing and 
back-scattered outgoing radiation, they find for the Reissner-Nordstrom (RN) 
solution that the mass function (qualitatively Rap-yS oc M/r 3 ) d iverg es at the 
CH (mass inflation). However, Ori showed both for RN and KerroEll that the 
metric perturbations are finite (even though K tivpa W lvt "' diverges) so that an 
observer would not be destroyed by tidal forces (the tidal distortion would be 
finite) and could survive passage through the CH. A numerical solution of the 
Einstein-Maxwell-scalariield equations could test these perturbative results. 

Gnedin and Gnedin E3 have numerically evolved the spherically symmetric 
Einstein-Maxwell with massless scalar field equations in a 2 + 2 formulation. 
The initial conditions place a scalar field on part of the RN event horizon 
(with zero field on the rest). An asymptotically null or spacelike singularity 
whose shape depends on the strength of the initial perturbation replaces the 
CH. For a sufficiently strong perturbation, the singularity is Schwarzschild- 
like. Although they claim to have found that the CH evolved to become a 
spacelike singularity, the diagrams in their paper show at least part of the final 
singularity to be null or asymptotically null in most cases. 

Brady and Smith H used the Goldwirth-Piran formulationS to study the 
same problem. They assume the spacetime is RN for v < vq. They follow 
the evolution of the CH into a null singularity, demonstrate mass inflation, 
and support (with observed exponential, decay of the metric component g) 
the validity of previous analytic results c3u3'E3Ea including the "weak" nature 
of the singularity that forms. They find that the observer hits the null CH 
singularity before falling into the curvature singularity ai_r = 0. Whether or 
not these results are in conflict with Gnedin and Gnedinl£3 is unclearE-H 

We finally mention Hiibner'sE3 numerical scheme to evolve on a conformal 
compactified grid using Friedrich's formalismEa He considers the spherically 
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symmetric scalar field model in a 2 + 2 formulation. So far, this, code has been 
used to locate singularities and to identify Choptuik's scalingt3 



3 Singularities in cosmological models 

3.1 Mixmaster dynamics 

Belinskii, Khalatnikov, and LifshitzB (BKL) described the singularity approach 
of vacuum Bianchi IX cosmologies as an infinite sequence of Kasncr epochs 
whose indices change when the scalar curvature terms in Einstein's equations 
become important. They were able to describe the dynamics approximately 
by ajroap evolving a discrete set of parameters from one Kasner epoch to the 
nextQ'LD For example, the Kasner indices for the power law dependence of the 
anisotropic scale factors can be parametrized by a single variable u > 1. BKL 
determined that 

Un+1 = Uu n -l)- 1 ' 2 l< n u n <2 ■ (6) 

The subtraction in the denominator for 1 < u n < 2 yields the sensitivity to ini- 
tial conditions associated with Mixmaster dynamics (MD). MisnerH described 
the same behavior in terms of the model's volume and anisotropic shears. A 
multiple of the scalar curvature acts as an outward moving potential in the 
anisotropy plane. Kasner epochs become straight line trajectories moving out 
a potential corner while bouncing from one side to the other. A change of 
corner ends a BKL era when u — > (u — Numerical evolution of Einstein's 

equations was used to explore the accuracy of the BKL map as a descriptor of 
the dynamics as well as the implications of the maplIaMlIZI 

Later, the [BKL sensitivity to initial conditions was discussed in the lan- 
guage of chaosE§t2l However, the chaotic nature of Mixmaster dynamics was 
questioned when numerical ciznliitirm of the Mixmaster equations yielded zero 
Lyapunov exponents (LE's)HE3£j (The LE measures the divergence of ini- 
tially nearby trajectories. Only an exponential divergence, characteristic of a 
chaotic system, will yield positive exponent.) Other numerical studies yielded 
positive LEEH This issue was resolved when the LE w as sho wn numerically and 
analytically to depend on the choice of time variableEH Although MD itse 



is well-understood, its c 
Recently, Le Blanc et al 



iracterization as chaotic or not is still controversial! 
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VEi have shown (analytically and numerically) that MD 
can arise in Bianchi VIo models with magnetic fields. In essence, the magnetic 
field provides the_wall needed to close the potential in a way that yields the 
BKL map for u. E-3 
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There are some recent numerical studies of Mixmasterjjiynamics in other 
theories of gravity. For example, Carretero-Gonzalez et alE3 find evidence p£ 
chaotic behavior in Bianchi IX-Brans-Dicke solutions while Cotsakis et alE3 
have shown that Bianchi IX models in 4th order gravity theories have stable 
non-chaotic solutions. 

3.2 Inhomogeneous Cosmologies 

In the remainder of this chapter, I shall discuss results of a joint project with V. 
Moncrief, D. GarfmklG,.aiid B. Grubisic to explore the nature of the generic cos- 
mological singularityE3E3 BKL have conjectured that one should expect such 
a singularity to be locally of the Mixmaster typeQ The main difficulty with 
the acceptance of this conjecture has been the ccptroversy payer whether the 
required time slicing pan be constructed globallyE3 Montanipl Belinskii(I2l and 
Kochnev and KirillovO have pointed out that if the BKL conjecture is correct, 
the spatial structure of the singularity could become extremely complicated as 
bounces occur at different locations at different times. The simplest cosmolog- 
ical models which might have local MD are vacuum universes on T 3 x R with a 
U(l) symmetryEj The non-commuting Killing vectors of local MD can be con- 
structed since only one Killing vector is already present. The two commuting 
Killing vectors of the even simpler plane symmetric Gowdy cosmologies BO 
preclude their use to test the conjecture. However, these models are interest- 
ing in their_pwn right since they have been conjectured to possess an AVTD 
singularity^ 

We employ a symplectic partial differential equation solverfaO a type of 
operator splitting which singles out the AVTD limit if it is present. If the 
model is, in fact, AVTD, the approximation in the numerical scheme should 
become more accurate as the singularity is approached. An outline of the 
method follows. 

For a field q(x, t) and its conjugate momentum n(x, t) split the Hamiltonian 
operator into kinetic and potential energy subhamiltonians. Thus, 

H = J dx{\v 2 + V[q]} =H 1 (TT)+H 2 (q). (7) 

If the vector X = (jr, q) defines the variables at time t, then the time evolution 
is given by 

dX , , . 

— = {H,X} PB =AX (8) 

where { }ps is the Poisson bracket. The usual exponentiation yields an 
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evolution operator 

e AAt = e A x (At/2) e A 2 At eAl (At/2) + q^S) (g) 

for A = A\ + A2 the generator of the time evolution. Higher order accuracy 
may be obtained by a better approximation to the evolution operatorlI3 This 
method is useful when exact solutions for the subhamiltonians are known. For 
the given H, variation of Hi yields the solution 

q = q + 7r At , 7T = 7r , (10) 

while that of H2 yields 

8V 

q = q , 7T = 7T - — At (11) 
Sq 

where SV /Sq is the appropriate functional derivative. Note that H 2 is exactly 
solvable for any potential V no matter how complicated although the required 
differenced form of the potential gradient may be non-trivial. One evolves from 
t to t + At using the exact solutions to the subhamiltonians according to the 
prescription given by the approximate evolution operator (||). 

The Gowdy model on T 3 x R sepjes as an excellent laboratory for this 
method. It is described by the metricO 

ds 2 = e- x / 2 e T / 2 (-e- 2T dT 2 +dd 2 ) 

+e- T [e p dcr 2 + 2e p Q da dS + {e p Q 2 + e~ p ) dS 2 } (12) 

where A, P, Q are functions of 9, r. We impose T 3 spatial topology by requiring 
< 0, a, 5 < 2ir and the metric functions to be periodic in 9. If we assume P 
and Q to be small, we find them to be respectively the amplitudes of the + and 
x polarizations of the gravitational waves with A describing the background in 
which they propagate. The time variable r measures the area in the symmetry 
plane with t = 00 a curvature singularity. Einstein's equations split into two 
groups. The first is nonlinear ly coupled wave equations for P and Q (where 
, a = d/da): 

P, TT -e- 2T P, ee -e 2P (Q 2 -e~ 2T Q 2 ) = 0, (13) 
Q, TT -e- 2T Q l8S + 2 (P, r Q, T -e- 2T P, e Q, e ) = 0. (14) 

The second contains the Hamiltonian and ^-momentum constraints respec- 
tively which can be expressed as first order equations for A in terms of P and 
Q: 

A JT - [P 2 + e- 2T P 2 + e 2P {Q 2 + e~ 2T Q 2 )] = 0, (15) 
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X,g- 2{P, g P, T + e 2P Q, e Q, T ) = 0. 



(16) 



This break into dynamical and constraint equations removes two of the most 
problematical areas of numerical relativity from this model. (1) The normally 
difficult initial value problem becomes trivial since P, Q and their first time 
derivatives may be specified arbitrarily (as long as the total 9 momentum in the 
waves vanishes). (2) The constraints, while guaranteed to be preserved in an 
analytic evolution by the Bianchi identities, are not automatically preserved in 
a numerical evolution with Einstein's equations in differenced form. However, 
in the Gowdy model, the constraints are trivial since A may be constructed 
from the numerically determined P and Q. For the special case of the polarized 
Gowdy modeL-(/2 = 0), P satisfies a linear wave equation whose exact solution 
is well- known E.3 For this case, it has been proven that the singularity is_AVTDEl 
This has also been conjectured to be true for generic Gowdy modelsLJ 

The wave equations (|l3|) can be obtained by variation of the Hamiltonian 




so that implementation of the symplectic method is straightforward. Again 
equations obtained from the separate variations of Hi and H 2 are exactly 
solvable. Variation of Hi yields the terms in ( |l3| ) containing time derivatives. 
These have the exact (AVTD) solution 

P = -j3r + ln[a (1 + C 2 e 2/3r )] -» /3t as r -> oo, 

c, —> Qa as t — ► oo , 



a(l + C 2 e 2/3r ) 

g(i - ee 2fiT ) 

(1 + C 2 e 2 ^) 

n Q - -2a(3<: (18) 



-/3(1-C 2 e 2 ^) 

nP = (1 + C 2 e 2/3r) ^ P as T -> OO , 



in terms of four constants a, /3, Q, and £. To employ the AVTD solution in 
the symplectic method, the values of P, Q, irp, and 7Tq at t j are used to find 
a, (3, (, and £. These are substituted in JIs| ) to evolve to new values at r J+1 
according to (||). Evolution with Hi is still easy because P and Q are constant. 
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For completeness, we give the (2nd order) differenced form of H2 as 



-2t 



Ho 



N 



(Qi — Qi-i) 



(19) 



The exponential prefactor e _2r in H 2 makes plausible the conjectured AVTD 
singularity. However, P — > vt (for v > 0) as r — > 00, (where from ( |l8| ) v = (3). 
If v > 1, the term e~ 2r e 2P Q,g in ( |l7[ ) can grow rather than decay as r — > 00. 
This has led to the conjecture that the .AVTD limit requires v < 1 everywhere 
except, perhaps, at isolated—values of #L3 

Our numerical results E3'E3 use the (standard) initial data P = 0, up = 
Vq cos 6, Q = cos 8, and ttq = 0. This model is actually generic for the following 
reasons: The cos (9 dependence is the smoothest nontrivial possibility. With 
cos n6, the solution is repeated n times on the grid yielding the same result 
with poorer resolution. The amplitude of Q is irrelevant since the Hamiltonian 
( |l7| ) is invariant under Q — > pQ, P —> P — hip. This also means that any 
unpolarized model is qualitatively different from a polarized (Q = 0) one no 
matter how small Q is. 

The accuracy and stability of the code easily allow verification of the con- 
jectured AVTD behaviorEll A plot of the maximum value of v vs r (Fig. [[]) 
shows strong support for the conjecture that v < 1 in the AVTD regime. 
However, Fig. [l] also shows that a simulation at higher spatial resolution be- 
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Figure 1: Plot of v max vs t. The maximum value oiv is found for two simulations with 3200 
(solid line) and 20000 spatial grid points (circles) respectively. The horizontal line indicates 
v = 1. 



gins to diverge from one at lower resolution. Normally, failure of convergence 
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signals numerical problems. Here something different is occurring. The evo- 
lution of spatial structure in P depends on competition between the two non- 
linear terms in the P wave equation. Approximating the wave equation by 
P, TT + either of the nonlinear terms = yields a first integral. The two po- 
tentials are V\ = 7Tge _2P and V% — e~ 2r e 2P Q,g. Non-generic behavior can 
arise at isolated points where either Q,g or ttq vanishes. Say such a point is 9q. 
The finer the grid, the closer will be some grid point to 9q. Thus non-generic 
behavior will become more visible on a finer grid. Detailed examination shows 
that the differences seen in Fig. [j] are due to the slower decay of v to a value 
below unity at isolated grid points in the higher resolution simulation. 

Details of the high resolution simulation are shown in Fig. p[ The narrow 




Figure 2: P (dashed line) and Q (solid line) vs 9 at r = 12.4 for the standard initial data set 
with do = 5 for < < 7r for a simulation containing 20000 spatial grid points in the interval 
[0, 2-7r]. The peaks in P are essentially the same in that they occur where Q,g ps while the 
apparent discontinuity in Q occurs where ttq RJ and P < (P = is the horizontal line). 

peaks in P occur where Q,e~ 0. Generically, if P f=s vt and v > 1, the 
potential V2 dominates. The relevant first integral of (|l^) is 

— J +e 2Z Q, 2 g = const (20) 

where Z = P — t. A bounce off V2 yields dZ/dr — > —dZ/dr or v — 1 — » 1 — v. 
If the new v < 0, then V\ dominates yielding the first integral 

P 2 +e- 2P TT 2 Q = const (21) 

causing P lT — * —P, T orm —v. Eventually, bouncing between potentials gives 
V < 1. However, if Q,g w 0, but is not precisely zero, it takes a long time for 
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the bounce off V2 to occur. Precisely at #0 (where Q,e~ 0), v > 1 persists. 
The apparent discontinuity in Q in Fig. ^|is not a numerical artifact. It occurs 
where P < and ttq w 0. Since Q, T — e~ 2P nQ and ttq w c(9 — 61) (if ttq = 
at 0i ), Q, r grows exponentially in opposite directions about 9\. The potential 
V\ drives P to positive values unless ttq = 0. Thus this feature will narrow as 
the simulation proceeds. 

Finally, we note that the spatial structure in P scales with the parameter 
vo in the initial data. In our standard initial data set, greater values of vq 
lead to the appearance of additional spatial structure in a shorter time. The 
rate of structure formation decreases and then stops as the AVTD regime is 
approached. The scaling is best for the time T5 at which the 5th peak appears 
in P although it is also seen for T3 and tj (the even nature of the solution 
causes two peaks to appear at once except at 9 = tt) and is described by 
T5 = 0(^0—^0°) _1 where, if vo = v^ 3 , the 5th peak does not appear until t = 00 
for some constant slope a. Explanation of this scaling is still in progress. 

Given the success of the symplectic method in studying the singularity be- 
havior of the Gowdy model, we can consider its extension to the case of U(l) 
symmetric cosmologies. Moncrief has shown 113 that cosmological models on 
T 3 x R with a spatial U(l) symmetry can be described by five degrees of free- 
dom {x, z, A, ip, uj} and their respective conjugate momenta {Pa;,p z ,PA,p, r}. 
All variables are functions of spatial variables u, v and time, r. If we define a 
conformal metric g a b in the u-v plane as g a b — c A e a b(x, z) where 



(22) 



has unit determinant and choose the conformal lapse N = e , Einstein's equa- 
tions can be obtained by variation of 

^dudv{±p 2 z + \e iz pl + ip 2 + ±e 4 V - \p\ - 2 PA 




_-2r r I _A _ab\ I _A _ab\ a , _A/_— 2z\ 

,b] 



+ e' 2T [(eV b ) , ab - (eV b ) , a A, b + e A (e~ 2z ) , [a x„ 



+ 2e A e Q V,a <p,b + ±e A e-^e ab u;, a w, 6 ] } 
= Hi+Hz, (23) 

Note particularly that 

Hi =H^(-2z,x) + H^(-2(p,uj) + H[(A) (24) 

where H^(P, Q) is the kinetic part of the Gowdy Hamiltonian (|l7|). Of course, 
Hf is just a free particle Hamiltonian for the degree of freedom associated 
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with A. This means that not only are the equations from H\ exactly solvable 
but also that the Gowdy coding can be used with essentially no change. The 
potential term Hi is very complicated. However, it still contains no momenta 
so its equations are trivially exactly solvable. Thus, at least in principle, the 
extension of the Gowdy code to the two spatial dimensions of the U(l) code is 
completely straightforward. 

There are three complications, however, which cause the U(l) problem 
to be more difficult. The first involves the initial value problem (IVP) — the 
constraints must be satisfied on the initial spacelike slice. The constraints are 

Ho = H - 2p A = (25) 

(where Ti. is the density in (p3|)) and 

H a = -2it a . b + paA j(I -pA, a +pip, a +ruj, a = (26) 

where 7r a is in the 2-space with metric e ab and is linear in p x and p z with each 
term containing one or the other. Moncrief has proposed a particular solution 
to the IVP. First, identically satisfy H a — by choosing 

Px = Pz — <p, a = w,o = ; pA = ce A (27) 

for c a constant. Then, solve Hq = for either r or p. Solution is possible 
for c > Cram such that r 2 or p 2 > 0. This allows x, z, A, and p or r to be 
freely specified. (Without loss of generality, it is possible to set x = z = 
initially to yield e a b flat. Such a condition may always be imposed at one time 
by rescaling u and v.) 

The second difficulty also involves the constraints. While the Bianchi iden- 
tities guarantee the preservation of the constraints by the Einstein evolution 
equations, there is no such guarantee for differenced evolution equations. At 
this stage of the project, we monitor the maximum value of the constraints 
vs r over the spatial grid but do nothing else to try to stay on the constraint 
hypersurface. 

The third difficulty, and the one that is proving to be the greatest obstacle, 
is instability associated with spatial differencing in two dimensions. In an 
attempt to control the instability, we have introduced a form of third order 
accurate spatial averaging. In both test cases and generic models, the averaging 
procedure has allowed the code to run longer. However, the averaging can lead 
to deviations of the numerical solution from the correct one. Fortunately, by 
comparing runs with and without averaging, these artifacts are easy to identify 
so that averaging can allow the code to run long enough for a conclusion about 
the asymptotic singularity behavior to be drawn. 
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Moncrief has provided a test case for the U (1) code. It again starts with a 
polarized Gowdy solution and transforms it as either a one-dimensional (9 — > u 
or v) or two-dimensional (6 — > f(u,v)) test problem to satisfy the U(l) equa- 
tions (including the constraints). As a one-dimensional example, the agreement 
is excellent and the code can be run to large r. Difficulties arise in the two- 
dimensional test problem in regions where the spatial derivatives are large. In 
application of the U(l) code to generic models, those which are AVTD can 
display nonlinear wave interactions before settling down to U — * (where U is 
the density corresponding to in ( p3| ) and thus contains all the terms with 
spatial derivatives), z,(p,A — > const r, and x,ui — > const. Increasing spatial 
resolution will yield narrower (and steeper) structures and thus may not help 
to cure instabilities due to steep gradients. 

Despite the limitations of the code, conclusions can be drawn for generic 
models in our restricted class of initial data. We shall consider the models as 
representative of subclasses of the data. Models with r = oj — are called 
polarized. This condition is compatible with the above solution to the IVP 
and is preserved identically by the (analytic and numerical) evolution equa- 
tions. Grubisic and Moncrief have conjectured that these polarized models are 
AVTDO Therefore, the first model is chosen to be polarized. It exhibits the 
conjectured AVTD behavior as shown in Fig. |^ for U. Detailed examination of 
the variables also shows AVTD behavior. Other polarized models also appear 
to be AVTD. Actually, this is not too surprising since a polarized model in our 
class of initial data must begin with p x — p z = r = lj = tp, a = so that only 
A, x, and z may be freely specified. But it is possible to set x = z — as well 
without loss of generality leaving only A arbitrary. Thus, polarized models 
in this class have a very simple structure. The second model is unpolarized 
with p given and the Hamiltonian constraint solved for r in the IVP. Spatial 
averaging increases stability and allows the model to be followed to the point 
where only artifacts have U ^ as is shown in Fig. [|. It certainly seems as 
if this model is also AVTD. Similar behavior is seen in other examples of this 
class. The last model has r given with p obtained by solving the Hamiltonian 
constraint in the IVP. Models of this type are less stable, probably due to the 
growth of a steep feature in u> which does not appear in the other cases. For 
this reason, the parameters must be kept small although the artifacts still ap- 
pear more prominently. Fig. ^ shows that U — > except near artifacts where 
no statement can be made. 

Thus, we may conclude that the application of the symplectic method to 
Einstein's equations for collapsing universes allows the nature of the singularity 
to be studied. Application to the Gowdy model has yielded strong support for 
its conjectured AVTD singularity and has allowed the discovery and study of 
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Figure 3: Frames of U(u,v,t) for the polarized model x = z = A = sinusini;, p\ = 
12e A , w = r = 0. Time increases to the right and downward. The final frame corresponds 
to U ~ everywhere. 

interesting small scale spatial structure and scaling. Further progress in under- 
standing the generic singularity of U(l) cosmologies requires imporvements in 
handling steep spatial gradients. (Spectral methods to evaluate spatial deriva- 
tives have been shown to work in the Gowdy simulations E3 and may yield 
improved behavior in the U(l) case.) Nevertheless there is strong support 
that (at least within our restricted class of initial data) polarized models are 
AVTD. There is also support for AVTD behavior in all generic models studied 
so far. Mixmaster-like bounces have not been seen while U appears to become 
small everywhere the results can be trusted. Several factors could account for 
this: (1) the BKL conjecture might be false; (2) the simulations have not run 
long enough; (3) Mixmaster behavior is present but hidden in our variables; or 
(4) our class of initial data is insufficiently generic. All these possibilities will 
be explored in studies in progress. 
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